%% Code for Section 6 of Acemoglu and Loebbing: Automation and Polarization
% This code replicates the robustness checks in Section 6.2, the
% counterfactuals in Section 6.3, the additional results in online
% Appendix B.5.6, and all figures and tables of the paper
% see online Appendix B.5 for an explanation of the numerical solution
% procedure of the model, the calibration, and the data sources

% specify path to base folder (Replication AL2025):
basefolder = fullfile('C:','Users','jonas.loebbing','Desktop','Replication AL2025');

% choose between full and basic execution:
% the basic execution loads the simulation results from the output folder
% and produces all tables and figures
% the full execution produces the simulation results in the output folder
% and creates all tables and figures
option = "full"; % set to "basic" or "full"

% choose whether figures shall be produced as tikzpictures or not
% this requires the matlab2tikz function available here:
% https://de.mathworks.com/matlabcentral/fileexchange/22022-matlab2tikz-matlab2tikz
optiontikz = "no"; % set to "yes" or "no"

% set default line width for figures
set(groot,'defaultLineLineWidth',2.0);

%% create parameter structure
p = parameters();

%% set external parameters
% elasticity of substitution between tasks, from Humlum (2021)
p.lambda = 0.5;

%% load calibration results
load(fullfile(basefolder,'Output','Calibration Results','calibration.mat')); % calibration without minimum wage

% set calibrated parameters
p.m = parammin(1);
p.n = parammin(2);
p.a = parammin(3);
p.b = parammin(4);
p.ak = parammin(5);
p.A = parammin(7);

% obtain q from sq (=parammin(6)) via q = min(qmax,qzero + sq*(qmax-qzero)), qmax s.t. q_2016<q_inf
qzero = targetmodelmin.details1980.qzero; % automation threshold
qinfinity = targetmodelmin.details1980.qinfinity; % infinite output threshold
qmax1980 = qinfinity + log(0.21)-2*p.qtol; % maximal q in 1980 s.t. q in 2016 is below q_infinity-p.qtol
q1980 = min(qmax1980,qzero + parammin(6)*(qmax1980-qzero)); % log capital productivity in 1980
q2016 = q1980 - log(0.21); % log capital productivity in 2016

% compute wage percentiles, log wages and log wage changes in 1980
% according to model
logw1980 = log(targetmodelmin.details1980.w); %log wages in 1980
sgrid = targetmodelmin.details1980.sgrid;
wperc1980 = (sgrid-targetmodelmin.details1980.sbar)./(1-targetmodelmin.details1980.sbar); %wage percentiles in 1980
logw2016 = log(interp1(...
    [targetmodelmin.details2016.sbottom(1:end-1);targetmodelmin.details2016.stop],...
    [targetmodelmin.details2016.wbottom(1:end-1);targetmodelmin.details2016.wtop],sgrid)); %log wages in 2016
deltalogw = logw2016-logw1980; %log wage changes
ind50 = find(sgrid>=0.5,1); %index of median wage
deltalogw50 = deltalogw(ind50); %median log wage change
ddeltalogw = deltalogw - deltalogw50; %log wage change relative to median wage change

%% load data
% note: all wage levels are measured in 2008 USD
% wage data from Output/Wage Data/wage.csv
% capital income share from Data/Capital Shares and User Costs/Acemoglu
% Loebbing Capital Shares and User Cost.xlsx
% minimum wage data from Data/Minimum Wage/Acemoglu Loebbing Minimum
% Wages.xlsx

% load wage data
wagedata = readtable(fullfile(basefolder,'Output','Wage Data','wage.csv'));

% extract variables from wage data
sdata = table2array(wagedata(:,1))./100; % wage percentile
logwdata = table2array(wagedata(:,2)); % log wage
deltalogwdata = table2array(wagedata(:,3)); % change in log wage due to automation, 1980-2016

% construct data moments used in calibration
% find indices of 10, 30, 50, 90 percentiles
ind10 = find(sdata>=0.1,1);
ind30 = find(sdata>=0.3,1);
ind50 = find(sdata>=0.5,1);
ind90 = find(sdata>=0.9,1);

% log wage changes due to automation relative to median log wage change
deltalogwdata50 = deltalogwdata(ind50);
ddeltalogwdata = deltalogwdata-deltalogwdata50;

% construct moments
logw5010data = logwdata(ind50)-logwdata(ind10);
logw9050data = logwdata(ind90)-logwdata(ind50);
alphakdata = 0.16; % income share of non-residential equipment and software, excl. all other capital
deltalogw3010data = deltalogwdata(ind30)-deltalogwdata(ind10);
deltalogw5030data = -ddeltalogwdata(ind30);
deltalogw9050data = ddeltalogwdata(ind90);
deltaalphakdata = 0.01; % change in income share of non-residential equipment and software, excl. other capital, 1980-2016
logw50data = logwdata(ind50);
% collect moments
data = [logw5010data; logw9050data; alphakdata; deltalogw3010data; deltalogw5030data; deltalogw9050data; deltaalphakdata; logw50data];

% minimum wage levels (in 2008 USD)
mw = 7.1; % effective minimum wage 1980
fedmw = 7.1; % federal minimum wage 1980
mw16 = 7.46; % effective minimum wage 2016
fedmw16 = 6.57; % federal minimum wage 2016
mw22 = 7.96; % effective minimum wage 2022
fedmw22 = 5.58; % federal minimum wage 2022 (7.25 USD in 2022)
maxmw22 = 12.38; % max state minimum wage 2022 (16.1 USD in 2022)

%% create tables and figures with calibration results (Section 6.1)
% Table 1A: parameter table
parameter = ["lambda"; "A"; "m"; "n"; "a"; "ak"; "b"; "log(q)"];
value = [p.lambda; p.A; p.m; p.n; p.a; p.ak; p.b; q1980];
parametertab = table(parameter,value);
writetable(parametertab, fullfile(basefolder,'Output','Tables','Table1a.txt'), 'Delimiter', '\t');

% Table 1B: comparison of data and model moments
moment = ["log(w_0.5)-log(w_0.1), 1980"; "log(w_0.9)-log(w_0.5), 1980";...
    "alpha_k"; "Delta (log(w_0.3)-log(w_0.1))";...
    "Delta (log(w_0.5)-log(w_0.3))"; "Delta (log(w_0.9)-log(w_0.5))";...
    "Delta alpha_k"; "log(w_0.5), 1980"];
model = targetmodelmin.targets(1:numel(data));
momenttab = table(moment,data,model);
writetable(momenttab,fullfile(basefolder,'Output','Tables','Table1b.txt'),'Delimiter','\t');

% Figure 4: log wage distribution in 1980, data vs model
% smooth data
logwdatasmooth = smoothdata(logwdata,"Gaussian",15);
% create figure
figlogwage80 = figure('Units','centimeters','Position',[2 1 16 9]);
plot(sdata,logwdatasmooth,'b','DisplayName','data');
hold on
plot(wperc1980,logw1980,'--r','DisplayName','model');
xlim([0 1]);
xlabel('Wage Percentile');
ylabel('Log Wage');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figlogwage80,fullfile(basefolder,'Output','Figures','fig4-logwage80.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig4-logwage80.tex'));
end

% Figure 5a: change in log wages due to automation, 1980-2016
figdeltalogwage = figure('Units','centimeters','Position',[2 1 9 9]);
yyaxis left
plot(sdata,deltalogwdata,'b','DisplayName','data (left axis)');
ylim([-0.1 0.25]);
yticks(linspace(min(ylim), 0.2, 4));
ylabel('Log Wage Change');
hold on
yyaxis right
plot(wperc1980,deltalogw,'--r','DisplayName','model (right axis)');
xlim([0 1]);
ylim([0.2 0.55]);
yticks(linspace(min(ylim), 0.5, 4));
xlabel('Wage Percentile');
legend('Location','northwest');
legend('boxoff');
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(1).Label.Color = 'k';
ax.YAxis(2).Color = 'r';
exportgraphics(figdeltalogwage,fullfile(basefolder,'Output','Figures','fig5a-deltalogwage.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig5a-deltalogwage.tex'));
end

% Figure 5b: change in assignment due to automation, 1980-2016
figassignment = figure('Units','centimeters','Position',[2 1 9 9]);
plot(targetmodelmin.details1980.x,targetmodelmin.details1980.sgrid,...
    'b','DisplayName','1980 capital cost');
hold on
plot(targetmodelmin.details2016.x,targetmodelmin.details2016.sgrid,...
    '--r','DisplayName','2016 capital cost');
xlim([0 1]);
ylim([0 1]);
xlabel('Task x');
ylabel('Wage Percentile');
legend('Location','southeast');
legend('boxoff');
exportgraphics(figassignment,fullfile(basefolder,'Output','Figures','fig5b-assignment.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig5b-assignment.tex'));
end

%% robustness check: elasticity of substitution between tasks (Section 6.2)

% create vector of alternative values for lambda
lambdavec = [0.3; 0.7];

% compute calibration targets for alternative values of lambda
if option == "full"

    parfor i = 1:numel(lambdavec)
    
        warning('off','all');
    
        % set lambda
        parrob = p;
        parrob.lambda = lambdavec(i);

        % compute targets
        targetmodelrob(i) = targetfun(parammin,0,0,parrob);

    end

    save(fullfile(basefolder,'Output','Simulation Results','robustness-lambda.mat'),'targetmodelrob','lambdavec');

else

    load(fullfile(basefolder,'Output','Simulation Results','robustness-lambda.mat'));

end

% Table 2, columns 1-4: robustness checks wrt lambda
modelrobustness = zeros(numel(data),numel(lambdavec)+1);
modelrobustness(:,1) = model;
for i = 1:numel(lambdavec)
    modelrobustness(:,i+1) = targetmodelrob(i).targets(1:numel(data));
end
momenttabrob = table(moment,data,modelrobustness);
writetable(momenttabrob, fullfile(basefolder,'Output','Tables','Table2(col1-4).txt'), 'Delimiter', '\t');

% compute 1980 log wages and 1980-2016 log wage change due to automation
% for different lambdas
% pre-allocation:
logw1980rob = cell(1,numel(lambdavec));
logw2016rob = cell(1,numel(lambdavec));
deltalogwrob = cell(1,numel(lambdavec));
% compute wages:
for i = 1:numel(lambdavec)
    % 1980 log wages:
    logw1980rob{i} = log(targetmodelrob(i).details1980.w);
    % 2016 log wages on 1980 worker grid:
    logw2016rob{i} = log(interp1([targetmodelrob(i).details2016.sbottom(1:end-1);...
        targetmodelrob(i).details2016.stop(1:end)],...
        [targetmodelrob(i).details2016.wbottom(1:end-1);...
        targetmodelrob(i).details2016.wtop(1:end)],...
        targetmodelrob(i).details1980.sgrid));
    % log wage change:
    deltalogwrob{i} = logw2016rob{i} - logw1980rob{i};
end

% Figure 6a: change in log wages due to automation, robustness check wrt
% lambda
figdeltalogwagerob = figure('Units','centimeters','Position',[2 1 9 9]);
yyaxis left
plot(sdata,deltalogwdata,'b','DisplayName','data (left axis)');
ylim([-0.15 0.4]);
ylabel('Log Wage Change');
hold on
yyaxis right
plot(wperc1980,deltalogw,'--r','DisplayName','lambda = 0.5 (right axis)');
plot(targetmodelrob(1).details1980.sgrid,deltalogwrob{1},'-.green',...
    'DisplayName','lambda = 0.3 (right axis)');
plot(targetmodelrob(2).details1980.sgrid,deltalogwrob{2},':black',...
    'DisplayName','lambda = 0.7 (right axis)');
xlim([0 1]);
ylim([0.15 0.7]);
xlabel('Wage Percentile');
legend('Location','northwest');
legend('boxoff');
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(1).Label.Color = 'k';
ax.YAxis(2).Color = 'r';
exportgraphics(figdeltalogwagerob,fullfile(basefolder,'Output','Figures','fig6a-deltalogwagerob.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig6a-deltalogwagerob.tex'));
end

%% robustness check: capital price index (Section 6.2)

% alternative capital price indices:
% eqqa: non-residential equipment, quality-adjusted
% sw: software
% baseline: non-residential equipment and software, quality-adjusted
% see "Acemoglu Loebbing Capital Shares and User Cost.xlsx" for information
% about price indices
qindexvec = ["eqqa"; "sw"];

% compute calibration targets using different price indices for the capital
% price decline between 1980 and 2016:
if option == "full"

    parfor i = 1:numel(qindexvec)
    
        warning('off','all');
    
        % choose price index:
        parrob = p;
        parrob.qindex = qindexvec(i);

        % compute targets
        targetmodelrobq(i) = targetfun(parammin,0,0,parrob);

    end

    save(fullfile(basefolder,'Output','Simulation Results','robustness-q.mat'),'targetmodelrobq','qindexvec');

else

    load(fullfile(basefolder,'Output','Simulation Results','robustness-q.mat'));

end

% Table 2, colums 1-2,5-6: robustness checks wrt capital price index
modelrobustnessq = zeros(numel(data),numel(qindexvec)+1);
modelrobustnessq(:,1) = model;
for i = 1:numel(qindexvec)
    modelrobustnessq(:,i+1) = targetmodelrobq(i).targets(1:numel(data));
end
momenttabrobq = table(moment,data,modelrobustnessq);
writetable(momenttabrobq, fullfile(basefolder,'Output','Tables','Table2(col1-2+5-6).txt'), 'Delimiter', '\t');

% compute 1980-2016 log wage change for different capital price declines
% pre-allocation:
logw1980robq = cell(1,numel(qindexvec));
logw2016robq = cell(1,numel(qindexvec));
deltalogwrobq = cell(1,numel(qindexvec));
% compute wages:
for i = 1:numel(qindexvec)
    % 1980 log wages:
    logw1980robq{i} = log([targetmodelrobq(i).details1980.wbottom(1:end-1);...
        targetmodelrobq(i).details1980.wtop]);
    % 2016 log wages on 1980 worker grid:
    logw2016robq{i} = log(interp1([targetmodelrobq(i).details2016.sbottom(1:end-1);...
        targetmodelrobq(i).details2016.stop(1:end)],...
        [targetmodelrobq(i).details2016.wbottom(1:end-1);...
        targetmodelrobq(i).details2016.wtop(1:end)],...
        [targetmodelrobq(i).details1980.sbottom(1:end-1);...
        targetmodelrobq(i).details1980.stop]));
    % log wage change:
    deltalogwrobq{i} = logw2016robq{i} - logw1980robq{i};
end

% Figure 6b: change in log wages due to automation, robustness check wrt
% capital price index
figdeltalogwagerobq = figure('Units','centimeters','Position',[2 1 9 9]);
yyaxis left
plot(sdata,deltalogwdata,'b','DisplayName','data (left axis)');
ylim([-0.15 0.4]);
ylabel('Log Wage Change');
hold on
yyaxis right
plot(wperc1980,deltalogw,'--r','DisplayName','baseline (right axis)');
plot([targetmodelrobq(1).details1980.sbottom(1:end-1);...
    targetmodelrobq(1).details1980.stop],deltalogwrobq{1},'-.green',...
    'DisplayName','equipment (right axis)');
plot([targetmodelrobq(end).details1980.sbottom(1:end-1);...
    targetmodelrobq(end).details1980.stop],deltalogwrobq{end},':black',...
    'DisplayName','software (right axis)');
xlim([0 1]);
ylim([0.15 0.7]);
xlabel('Wage Percentile');
legend('Location','northwest');
legend('boxoff');
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(1).Label.Color = 'k';
ax.YAxis(2).Color = 'r';
exportgraphics(figdeltalogwagerobq,fullfile(basefolder,'Output','Figures','fig6b-deltalogwagerobq.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig6b-deltalogwagerobq.tex'));
end

%% counterfactual 1: increasing capital productivity (Section 6.3)

% create vector of capital productivity levels, with steps equal to the
% increase in capital productivity from 1980 to 2016
qvec = q1980:-log(0.21):1;

% compute equilibrium for each productivity level
if option == "full"

    % use 1980 and 2016 equilibrium for first entries
    eqq(1) = targetmodelmin.details1980;
    eqq(2) = targetmodelmin.details2016;

    parfor i = 3:numel(qvec)
    
        warning('off','all')
    
        eqq(i) = solvemodelmw(p,qvec(i),0);
    
    end
    
    save(fullfile(basefolder,'Output','Simulation Results','counterfactual-capprod.mat'),"eqq","qvec");

else 

    load(fullfile(basefolder,'Output','Simulation Results','counterfactual-capprod.mat'));

end

% compute log wage changes for first counterfactual (starting from its 2016
% level, capital productivity increases by the same factor as from 1980 to
% 2016)
logwcf = log(interp1([eqq(3).sbottom(1:end-1);eqq(3).stop],...
    [eqq(3).wbottom(1:end-1);eqq(3).wtop],sgrid)); %counterfactual log wages
deltalogwcf = logwcf-logw2016; %log wage changes relative to 2016
deltalogwcf80 = logwcf-logw1980; %log wage changes relative to 1980

% Figure 7a: change in wages due to a counterfactual fall in capital costs,
% log wage changes
figdeltalogwagecfq = figure('Units','centimeters','Position',[2 1 9 9]);
plot(wperc1980,deltalogw,'b','DisplayName','fall in capital cost 1980-2016');
hold on;
plot(wperc1980,deltalogwcf,'--r','DisplayName','counterfactual');
xlim([0 1]);
ylim([0.15 0.55]);
yticks(linspace(0.2, 0.5, 4));
xlabel('Wage Percentile');
ylabel('Log Wage Change');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figdeltalogwagecfq,fullfile(basefolder,'Output','Figures','fig7a-deltalogwagecfq.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig7a-deltalogwagecfq.tex'));
end

% Figure 7b: change in wages due to a counterfactual fall in capital costs,
% wage distributions
figlogwagecfq = figure('Units','centimeters','Position',[2 1 9 9]);
plot(wperc1980,logw1980,':black','DisplayName','1980 capital cost');
hold on;
plot(sgrid,logw2016,'b','DisplayName','2016 capital cost');
plot(sgrid,logwcf,'--red','DisplayName','counterfactual capital cost');
xlim([0 1]);
ylim([1.5 5.5]);
yticks(linspace(2, 5, 4));
xlabel('Wage Percentile');
ylabel('Log Wage');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figlogwagecfq,fullfile(basefolder,'Output','Figures','fig7b-logwagecfq.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig7b-logwagecfq.tex'));
end

%% Counterfactual 2: artificial intelligence (Section 6.3)

% create vector of curvature levels for capital productivity
bvec = p.b:(-p.b/10):p.b*3/10;

% for each b, adjust ak and q such that the peak productivity of capital
% relative to the median worker remains unchanged: i.e., when reducing the
% curvature, the capabilities of capital become broader but not better at
% their peak
% location of peak productivity of capital relative to median worker:
%xhat = -p.ak/(2*p.b);
xhat = (0.5*p.a-p.ak)/(2*p.b);
% peak productivity:
pihat = -p.b*xhat^2+q2016;
% adjust ak:
akvec = p.a/2 - 2*bvec*xhat;
% adjust q:
qvecb = pihat + bvec*xhat^2;

% compute equilibrium for each parameter vector
if option == "full"
    
    % initial parameters:
    eqb(1) = eqq(2);

    % new parameter vectors:
    parfor i = 2:numel(bvec)
    
        warning('off','all');
    
        % replace parameters
        par(i) = p;
        par(i).b = bvec(i);
        par(i).ak = akvec(i);
        
        % compute equilibrium
        eqb(i) = solvemodelmw(par(i),qvecb(i),0);
    
    end

    save(fullfile(basefolder,'Output','Simulation Results','counterfactual-curvature.mat'),'eqb','bvec','akvec','qvecb');

else

    load(fullfile(basefolder,'Output','Simulation Results','counterfactual-curvature.mat'));

end

% compute log wage changes due to decrease in absolute value of b of 60%, which yields a
% similar net output gain as increasing q by the same factor as between
% 1980 and 2016
logwb = log(interp1([eqb(7).sbottom(1:end-1);eqb(7).stop],...
    [eqb(7).wbottom(1:end-1);eqb(7).wtop],sgrid)); %counterfactual log wages
deltalogwb = logwb-logw2016; %log wage changes relative to 2016

% compute log capital productivity relative to log median labor
% productivity
% adjust parameters
parb = p;
parb.b = bvec(7);
parb.ak = akvec(7);
parb.q = qvecb(7);
p.q = q2016; % include 2016 capital productivity in parameter structure
% compute change in capital productivity due to reduction in curvature
logpsib = logcapprod(0:0.01:1,parb)-logcapprod(0:0.01:1,p);

% Figure 8a: change in wages due to a counterfactual progress in artificial
% intelligence, log wage changes
figdeltalogwagecfb = figure('Units','centimeters','Position',[2 1 9 9]);
plot(wperc1980,deltalogwcf,'blue','DisplayName','uniform fall in capital cost');
hold on
plot(wperc1980,deltalogwb,'--red','DisplayName','artificial intelligence');
xlim([0 1]);
ylim([-0.2 0.7]);
yticks(linspace(-0.2, 0.6, 5));
xlabel('Wage Percentile');
ylabel('Log Wage Change');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figdeltalogwagecfb,fullfile(basefolder,'Output','Figures','fig8a-deltalogwagecfb.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig8a-deltalogwagecfb.tex'));
end

% Figure 8b: change in wages due to a counterfactual progress in artificial
% intelligence, wage distributions
figlogwagecfb = figure('Units','centimeters','Position',[2 1 9 9]);
plot(wperc1980,logw2016,':black','DisplayName','2016 capital cost');
hold on
plot(wperc1980,logwcf,'blue','DisplayName','uniform fall in capital cost');
plot(wperc1980,logwb,'--red','DisplayName','artificial intelligence');
xlim([0 1]);
ylim([1.5 5.5]);
yticks(linspace(2, 5, 4));
xlabel('Wage Percentile');
ylabel('Log Wage');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figlogwagecfb,fullfile(basefolder,'Output','Figures','fig8b-logwagecfb.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig8b-logwagecfb.tex'));
end

% Figure B-2: counterfactual capital productivity changes
figcapprodcfb = figure('Units','centimeters','Position',[2 1 16 9]);
plot(0:0.01:1,logpsib,'blue');
xlim([0 1]);
ylim([0 140]);
xlabel('Task x');
ylabel('Log change in capital productivity');
exportgraphics(figcapprodcfb,fullfile(basefolder,'Output','Figures','figb2-capprodcfb.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','figb2-capprodcfb.tex'));
end

%% counterfactuals 3 and 4: minimum wage increase and low-skill supply reduction (Section 6.3)

% minimum wage: going from 2022 effective minimum wage to the highest state-level minimum
% wage

% compute equilibrium for 2022 federal minimum wage, 2022 effective mw, and max state mw
if option == "full"

    warning('off','all');

    % federal minimum wage
    eqmwnew(1) = solvemodelmw(p,q2016,fedmw22);

    % effective minimum wage
    eqmwnew(2) = solvemodelmw(p,q2016,mw22);

    % max state-level minimum wage
    eqmwnew(3) = solvemodelmw(p,q2016,maxmw22);

    save(fullfile(basefolder,'Output','Simulation Results','counterfactual-mw-new.mat'),"eqmwnew","fedmw22","mw22","maxmw22");

else

    load(fullfile(basefolder,'Output','Simulation Results','counterfactual-mw-new.mat'));

end

% compute log wage changes
% log wage under federal mw on sgrid
logwfedmw = log(interp1([eqmwnew(1).sbottom(1:end-1);eqmwnew(1).stop],...
    [eqmwnew(1).wbottom(1:end-1);eqmwnew(1).wtop],sgrid));
% log wage under effective mw on sgrid
logwmw = log(interp1([eqmwnew(2).sbottom(1:end-1);eqmwnew(2).stop],...
    [eqmwnew(2).wbottom(1:end-1);eqmwnew(2).wtop],sgrid));
% log wage under max state mw on sgrid
logwmaxmw = log(interp1([eqmwnew(3).sbottom(1:end-1);eqmwnew(3).stop],...
    [eqmwnew(3).wbottom(1:end-1);eqmwnew(3).wtop],sgrid));
% log wage change effective to max
deltalogwmw = logwmaxmw - logwmw;
% log wage change federal to max
deltalogwmwfedmax = logwmaxmw - logwfedmw;

% compute wage percentiles with minimum wage
for i = 1:3
    eqmwnew(i).wperc = (eqmwnew(i).sgrid-eqmwnew(i).sbar)/(1-eqmwnew(i).sbar);
end

% reduction in low-skill labor supply

% indicate to use new labor supply vector given by function labornew.m
p.labor = "new";

% compute equilibrium under new labor supply
if option == "full"

    % initial labor supply:
    eql(1) = targetmodelmin.details2016;

    warning('off','all');

    % new labor supply:
    eql(2) = solvemodelmw(p,q2016,0);

    % compute labor supply vector on grid:
    laborvec = labornew(eql(2).sgrid);

    save(fullfile(basefolder,'Output','Simulation Results','counterfactual-laborsupply.mat'),'eql','laborvec');

else

    load(fullfile(basefolder,'Output','Simulation Results','counterfactual-laborsupply.mat'));

end

% compute change in labor supply in percent
deltalpercent = 100*(laborvec-ones(size(eql(2).sgrid)));

% compute log wage changes due to labor supply change
% map log wages after ls change initial worker grid
logwl = log(interp1([eql(2).sbottom(1:end-1);eql(2).stop],...
    [eql(2).wbottom(1:end-1);eql(2).wtop],sgrid));
deltalogwl = logwl-logw2016; %log wage changes

% compute wage percentiles after labor supply change
eql(2).wperc = cumtrapz(eql(2).sgrid,laborvec)/trapz(eql(2).sgrid,laborvec);

% reset labor supply to baseline
p.labor = "baseline";

% Figure 9a: change in wages due to minimum wage increase and low-skill
% supply reduction, log wage changes
figdeltalogwagecfmw = figure('Units','centimeters','Position',[2 1 9 9]);
plot(sgrid,deltalogwmw,'b','DisplayName','minimum wage increase');
hold on
plot(sgrid,deltalogwl,'--red','DisplayName','low-skill supply reduction');
xlim([0 1]);
ylim([-0.06 0.17]);
xlabel('Wage Percentile');
ylabel('Log Wage Change');
legend('Location','northeast');
legend('boxoff');
exportgraphics(figdeltalogwagecfmw,fullfile(basefolder,'Output','Figures','fig9a-deltalogwagecfmw.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig9a-deltalogwagecfmw.tex'));
end

% Figure 9b: change in wages due to minimum wage increase and low-skill
% supply reduction, wage distributions
figlogwagecfmw = figure('Units','centimeters','Position',[2 1 9 9]);
plot(sgrid,logw2016,':black','DisplayName','baseline');
hold on
plot(sgrid,logwmaxmw,'blue','DisplayName','max state minimum wage');
plot(sgrid,logwl,'--red','DisplayName','reduced low-skill supply');
xlim([0 1]);
ylim([1.5 5.5]);
yticks(linspace(2, 5, 4));
xlabel('Wage Percentile');
ylabel('Log Wage');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figlogwagecfmw,fullfile(basefolder,'Output','Figures','fig9b-logwagecfmw.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig9b-logwage-cfmw.tex'));
end

% Figure B-3: counterfactual labor supply change
figdeltalabor = figure('Units','centimeters','Position',[2 1 16 9]);
plot(eql(2).sgrid,deltalpercent,'blue','DisplayName','labor supply change');
xlim([0 1]);
ylim([-10 0]);
xlabel('Initial Wage Percentile');
ylabel('Labor Supply Change in %');
exportgraphics(figdeltalabor,fullfile(basefolder,'Output','Figures','figb3-deltalabor.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','figb3-deltalabor.tex'));
end

% Figure B4a: change in assignment due to minimum wage increase
figassignmentcfmw = figure('Units','centimeters','Position',[2 1 9 9]);
plot(eqmwnew(2).x,eqmwnew(2).wperc,...
    'blue','DisplayName','effective minimum wage');
hold on
plot(eqmwnew(3).x,eqmwnew(3).wperc,'--red','DisplayName','max state minimum wage');
xlim([0 1]);
ylim([0 1]);
xlabel('Task x');
ylabel('Wage Percentile');
legend('Location','southeast');
legend('boxoff');
exportgraphics(figassignmentcfmw,fullfile(basefolder,'Output','Figures','figb4a-assignmentcfmw.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','figb4a-assignmentcfmw.tex'));
end

%Figure B4b: change in assignment due to low-skill supply reduction
figassignmentcfl = figure('Units','centimeters','Position',[2 1 9 9]);
plot(eql(1).x,eql(1).sgrid,...
    'blue','DisplayName','baseline');
hold on;
plot(eql(2).x,eql(2).wperc,'--red','DisplayName',...
    'reduced low-skill supply');
xlim([0 1]);
ylim([0 1]);
xlabel('Task x');
ylabel('Wage Percentile');
legend('Location','southeast');
legend('boxoff');
exportgraphics(figassignmentcfl,fullfile(basefolder,'Output','Figures','figb4b-assignmentcfl.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','figb4b-assignmentcfl.tex'));
end

%% Calibration results with minimum wage (online appendix B.5.7)
% load calibration results
load(fullfile(basefolder,'Output','Calibration Results','calibrationmw.mat'));

% set calibrated parameters
p.m = paramminmw(1);
p.n = paramminmw(2);
p.a = paramminmw(3);
p.b = paramminmw(4);
p.ak = paramminmw(5);
p.A = paramminmw(7);

% obtain q from sq (=parammin(6)) via q = min(qmax,qzero + sq*(qmax-qzero)), qmax s.t. q_2016<q_inf
qzeromw = targetmodelminmw.details1980.qzero; % automation threshold
qinfinitymw = targetmodelminmw.details1980.qinfinity; % infinite output threshold
qmax1980mw = qinfinitymw + log(0.21)-2*p.qtol; % maximal q in 1980 s.t. q in 2016 is below q_infinity-p.qtol
q1980mw = min(qmax1980mw,qzeromw + paramminmw(6)*(qmax1980mw-qzeromw)); % log capital productivity in 1980
q2016mw = q1980mw - log(0.21); % log capital productivity in 2016

% compute wage percentiles, log wages and log wage changes in 1980
% according to model
logw1980mw = log(targetmodelminmw.details1980.w); %log wages in 1980
sgridmw = targetmodelminmw.details1980.sgrid;
wperc1980mw = (sgridmw-targetmodelminmw.details1980.sbar)./(1-targetmodelminmw.details1980.sbar); %wage percentiles in 1980
logw2016mw = log(interp1(...
    [targetmodelminmw.details2016.sbottom(1:end-1);targetmodelminmw.details2016.stop],...
    [targetmodelminmw.details2016.wbottom(1:end-1);targetmodelminmw.details2016.wtop],sgridmw)); %log wages in 2016
deltalogwmw = logw2016mw-logw1980mw; %log wage changes
ind50mw = find(wperc1980mw>=0.5,1); %index of median wage
deltalogw50mw = deltalogw(ind50mw); %median log wage change
ddeltalogwmw = deltalogwmw - deltalogw50mw; %log wage change relative to median wage change

% Table B-1A: Calibrated Parameters (with minimum wage)
parametermw = ["lambda"; "A"; "m"; "n"; "a"; "ak"; "b"; "log(q)"];
valuemw = [p.lambda; p.A; p.m; p.n; p.a; p.ak; p.b; q1980];
parametertabmw = table(parametermw,valuemw);
writetable(parametertabmw, fullfile(basefolder,'Output','Tables','Tableb1a.txt'), 'Delimiter', '\t');

% Table B-1B: Comparison of Moments (with minimum wage)
momentmw = ["log(w_0.5)-log(w_0.1), 1980"; "log(w_0.9)-log(w_0.5), 1980";...
    "alpha_k"; "Delta (log(w_0.3)-log(w_0.1))";...
    "Delta (log(w_0.5)-log(w_0.3))"; "Delta (log(w_0.9)-log(w_0.5))";...
    "Delta alpha_k"; "log(w_0.5), 1980"];
modelmw = targetmodelminmw.targets(1:numel(data));
momenttabmw = table(momentmw,data,modelmw);
writetable(momenttabmw, fullfile(basefolder,'Output','Tables','Tableb1b.txt'), 'Delimiter', '\t');

% Figure B-5a: Log wage distribution in 1980, data vs model (with minimum
% wage)
% smooth data
logwdatasmooth = smoothdata(logwdata,"Gaussian",15);
% create figure
figlogwage80mw = figure('Units','centimeters','Position',[2 1 9 9]);
plot(sdata,logwdatasmooth,'b','DisplayName','data');
hold on
plot(wperc1980mw,logw1980mw,'--red','DisplayName','model');
xlim([0 1]);
ylim([1.5 4]);
xlabel('Wage Percentile');
ylabel('Log Wage');
legend('Location','northwest');
legend('boxoff');
exportgraphics(figlogwage80mw,fullfile(basefolder,'Output','Figures','figb5a-logwage80mw.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','fig5ba-logwage80mw.tex'));
end

% Figure B-5b: change in log wages due to automation, data vs model (with
% minimum wage
figdeltalogwagemw = figure('Units','centimeters','Position',[2 1 9 9]);
yyaxis left
plot(sdata,deltalogwdata,'b','DisplayName','data (left axis)');
ylim([-0.1 0.25]);
yticks(linspace(min(ylim), 0.2, 4));
ylabel('Log Wage Change');
hold on
yyaxis right
plot(wperc1980mw,deltalogwmw,'--r','DisplayName','model (right axis)');
xlim([0 1]);
ylim([0.2 0.55]);
yticks(linspace(min(ylim), 0.5, 4));
xlabel('Wage Percentile');
legend('Location','northwest');
legend('boxoff');
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(1).Label.Color = 'k';
ax.YAxis(2).Color = 'r';
exportgraphics(figdeltalogwagemw,fullfile(basefolder,'Output','Figures','figb5b-deltalogwagemw.pdf'));
if optiontikz == "yes"
    matlab2tikz(fullfile(basefolder,'Output','Figures','figb5b-deltalogwagemw.tex'));
end